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ABSTRACT 

The early evolution of dense star clusters is possibly dominated by close interac- 
tions between stars, and physical collisions between stars may occur quite frequently. 
Simulating a stellar collision event can be an intensive numerical task, as detailed 
calculations of this process require hydrodynamic simulations in three dimensions. We 
present a computationally inexpensive method in which we approximate the merger 
process, including shock heating, hydrodynamic mixing and mass loss, with a simple 
algorithm based on conservation laws and a basic qualitative understanding of the 
hydrodynamics of stellar mergers. The algorithm relies on Archimedes' principle to 
dictate the distribution of the fluid in the stable equilibrium situation. We calibrate 
and apply the method to mergers of massive stars, as these are expected to occur in 
young and dense star clusters. We find that without the effects of microscopic mixing, 
the temperature and chemical composition profiles in a collision product can become 
double- valued functions of enclosed mass. Such an unphysical situation is mended by 
simulating microscopic mixing as a post-collision effect. In this way we find that head- 
on collisions between stars of the same spectral type result in substantial mixing, while 
mergers between stars of different spectral type, such as type B and O stars (~10 and 
^4OM0respectively), are subject to relatively little hydrodynamic mixing. 



1 INTRODUCTION 

The cores of star clusters can become so dense that close 
encounte rs and even collisions betwe en stars become quite 
common (|Portegies Zwart et al.|[l999l ). This effect is partic- 
ularly important in young star clusters where a broad range 
of stellar masses dominates the dynamical evolution of the 
cluster core. In a star cluster with an initial half-mass relax- 
ation time smaller than a few tens of million years, core col- 
lapse can occur well before the most massive stars leave the 
main-sequence and explode in a supernova. In such star clus- 
ters, these massive stars will dominate the core dynamics. 
This is most noticeable in both the arrest of core collapse by 
the formation of tight massive binaries and the occurrence 
of physical c ollisions between sta rs which are mediated by 
the binaries (jGaburov et al.|[2007h . 

Simulating stellar collisions in globular clusters has at- 
tracted considerable attention over the last decade. These 
calculations were mainly initiated to e xplain the pres- 
ence and characterist i cs of b lue stragglers IjSills et alj|2002t 
ILombardi et al. 1 12001 [ 19961 ). In many cases, collisions oc- 
cur infrequently enough that the already computationally 
expensive dynamical modelling of a cluster would hardly 
be slowed down by a simple treatment of collisions. How- 
ever, it was recently recognised that the number of col- 
lisions in young star clusters can become quite large, 
and this may lead to formation of a very mass ive object 
|Portegies Zwart et al. 1 120041 ; iFreitag et all 12006ft . Further- 
more, the subsequent evolution of such objects, and of colli- 



sion products in general, can depend on their internal struc- 
ture. 



To determine the structure of a collision product ac- 
curately, one can resort to high-resolution numerical sim- 
ulations of stellar mergers. Such simulations are generally 
carried out by mea ns of the smoothe d particle hydrodynam- 
ics (SPH) method (jMonaghan 2005). This however becomes 
computationally expensive, especially if one is interested in 
multiple collisions which may occur in the course of iV-body 
simulations. On the other hand, components of multiple sys- 
tems, such_as_Wnariesjuid triples, may experience a colli- 
sion IjGaburov et al.ll2007l ). and careful treatment of such 
collisions might play an important role for further evolu- 
tion of these multiple systems. Therefore, faster methods are 
needed to compute the outcome of a stellar merger event, 
and these methods can eventually be included in iV-body 
simulations. 



ILombardi et al l (|l996l ) constructed a method which has 
been successfully used for collisions between low-mass main- 
sequence stars, and in this letter we extend this method to 
collisions between massive stars. Our method has a general 
applicability, in the sense that it is suitable for mergers of 
all stellar types, including compact objects. 
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2 METHODS 

2.1 Guiding Principle 

The guiding idea of our method is Archimedes' principle. 
To determine how fluid will redistribute itself when brought 
out of equilibrium, perhaps by a stellar merger, we need 
to consider densities: a fluid element with a greater density 
than its environment will accelerate downward, while one 
with a smaller density will be buoyed upward. Any such 
element can ultimately settle into hydrostatic equilibrium 
once its density equals that of its environment. The density 
of a fluid element will not in general remain constant during 
the collision. Instead, it will be continuously adjusted in such 
a way that pressure equilibrium with the environment is 
achieved. In other words, parcels of gas expand or contract as 
necessary to equilibrate their pressure with the environment. 

We begin by considering a distribution of gas, in our 
case stars, and specify distribution functions for the pres- 
sure P, density p, and abundance of chemical species Xi 
as a function of the enclosed mass. For a given equation of 
state, the expression for the entropy can be determined. It is 
often mathematically simpler, although not formally neces- 
sary, to define an entropic variable A which is related to the 
specific entropy and conserved by fluid elements in adiabatic 
processes; we call this entropy variable A the buoyancy. 

The buoyancy can be calculated from the equation of 
state, P — P(A, p, Xi). In the case of an ideal gas A = 
P/p T , where F is the adiabatic index. By construction, the 
buoyancy A depends directly upon entropy and composition, 
and it remains constant for each fluid element in the absence 
of heating and mixing. Therefore, it is an important starting 
point for understanding the hydrodynamics of collisions. 

This idea has been successfully used to describe the un- 
derlying hydro dynamics of collisions among low mass main 
sequence stars (|Lombardi et al.lll99 6, 2002, 2003), which are 
well described by a monatomic ideal gas equation of state. 
In this case p = (P/A) 3 / 5 , such that at a given pressure, 
a smaller A yields a greater value of p. Consequently, fluid 
with smaller values of A sink to the bottom of a potential 
well, and the A profile of a star in stable hydrostatic equi- 
librium increases radially outwards. Indeed, it is straightfor- 
ward to show that the condition dA/dm > is equivalent 
to the usual Led oux criterion for conv ective stability of a 
nonrotating star (|Lombardi et al.lll996l ). 



2.2 Sorting Algorithm 

The fluid in high mass main-sequence stars is well described 
by a mixture of monatomic ideal gas and radiation in ther- 
mal equilibrium. In such cases, the total pressure is 



in Eq.[2j we define the buoyancy A through the following 
relation 
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where k is the Boltzmann constant, T is temperature, p is 
the mean molecular m ass, and a is th e radiation constant. 
The specific entropy is l|Clavtonll 19831 ) 
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where the quantity so depends only on composition and p 
is the ratio of gas to total pressure. To achieve the equality 
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Through manipulations with Eq. [T] Eq. [2] Eq. [3] and defini- 
tion of p, we find 
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Together, Eq. [3] and Eq. [4] provide the desired relationships 
between the pressure, buoyancy, density, and composition. 
Given a prescription for shock heating and an estimate of 
the mass loss, it is possible to determine the structure of 
a non-rotating product resulting from a head-on collision of 
two massive main-sequence stars, which we call parent stars. 

To begin with, we consider fluid elements at the centre 
of each parent star and determine which one will settle in 
the centre of the collision product. For this purpose, we both 
make an initial estimate of the central pressure of the colli- 
sion product and test which fluid element would have larger 
density when brought to this central pressure: a process that 
can be completed in two steps. First, as the composition and 
post-shock buoyancy of each fluid element can be considered 
known, Eq. [3] gives what would be the final (3 value of each 
element if it were to settle at the centre of the collision prod- 
uct. Second, these /3 values correspond to densities through 
Eq. [4] The fluid element with the greater density is the one 
that actually settles in the centre of the collision product. 
Next, we integrate the equations of hydrostatic equilibrium 
to determine the local pressure P throughout the product. 
This local pressure is used to determine which of subsequent 
fluid elements from the parent stars contribute to the struc- 
ture of the product. This procedure is repeated until the 
mass supply of the parent stars is depleted. As in the case of 
low-mass main-sequence stellar collisions, the central pres- 
sure is iteratively improved until the outer boundary condi- 
tion, namely that the surface pressure vanishes, is satisfied. 

As a consequence of this merging procedure, one must 
expect that fluid of a certain composition that has origi- 
nated in one parent star can be located arbitrarily close to 
fluid of a different composition from the other parent. These 
fluid elements must be in pressure equilibrium and also, by 
the condition of hydrostatic equilibrium, have the same den- 
sity. Therefore, in order to have adjacent fluid elements of 
different p, Eq. [3] and Eq. [4] imply that these elements must 
have different buoyancies A. We therefore should not expect 
that the A and p profiles in a nascent collision product will 
be single- valued, at least not before mixing has occurred on 
the microscopic level. 

If the parent stars have essentially the same homoge- 
neous composition as in the merger of two unevolved stars, 
then p will be uniform throughout the collision product, and 
for any given pressure, greater densities will correspond to 
smaller values of the buoyancy. In other words, the struc- 
ture of the collision product can then be simply determined 
by sorting the fluid in order of increasing A: the lowest A 
fluid is placed at the centre of the product and is surrounded 
by shells with increasingly higher values of A. The merging 
procedure also reduces to sorting by A in another special 
case, namely when radiation pressure is negligible (P —> 1), 
as Eq. [3]clearly reduces to the monatomic ideal gas result in 
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this limit. The merging procedure presented here is there- 
fore a n atural generalization of that presented for low mass 
stars in lLombardi etafl (|2002l . l2003h . 



2.3 Stability Criterion in High Mass Stars 

The unusual profiles that can be achieved in stellar merg- 
ers suggest that we take a careful look at the condition for 
convective stability: 



< 



(5) 



Here, the subscripts "env" and "el" refer to the environment 
and the fluid element respectively. Since pressure equilib- 
rium between the element and its immediate environment 
is established nearly instantaneously, the condition for con- 
vective stability can be written in the following way: 
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Using Eq. [3] and Eq. [4] we evaluate derivatives on both sides 
of Eq. |6] We then note that 
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and this allows us to write the stability condition in the 
following form 
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As the dlogp/dlogm is always negative throughout the 
star, we can rewrite the stability condition in a more conve- 
nient form: 

A 5 
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In a star in which Eq. [5]is satisfied, a perturbed element will 
experience restoring forces that cause it to oscillate around 
its equilibrium position. It can be seen that in the limit of 
ideal gas (Ti = 5/3), or when the composition is uniform, we 
recover the usual stability criteria dA/dm > 0. In the gen- 
eral case, however, A will typically increase outward, but can 
decrease in regions with an inverted composition gradient. 



3 VALIDATION 
3.1 Initial conditions 

In the previous section we have presented a sorting algo- 
rithm that generates a collision product in hydrostatic equi- 
librium. As the method deals with non-rotating products, 
we validate it by carrying out simulations of head-on colli- 
sions between main-sequence stars of different mass and age. 
The full set of simulations is presented in Tab. [T] 

Simulations of stellar collisions were carrie d out by 
mean s of a modified version of the GADGET2 code (|Springell 
2005). In particular, we included radiation pressure in the 
equation of state and the functionality which allows gener- 
ation of rela xed models of parent s tars in quasi-hydro static 
equilibrium l|Lombardi et al.ll2006l : iTurner et al.lll995h . 

We have prepared three-dimensional SPH stellar models 
from one-dimensional stellar models computed with the EZ 



Model 


Mi 


M 2 


Evolutionary State 


N 


T88, H88 


80 


8 


TAMS, HAMS 


880k 


T48, H48 


40 


8 


TAMS, HAMS 


480k 


T28, H28, Z28 


20 


8 


TAMS, HAMS, ZAMS 


280k 


T18, HIS, Z18 


10 


8 


TAMS, HAMS, ZAMS 


180k 



Table 1. The simulations carried out in this work. The model 
name of the simulations are given in the first column. The masses 
of the primary and the secondary stars are shown in the sec- 
ond and the third column respectively. The evolutionary state of 
the parent stars is given in the fourth column: TAMS, HAMS, 
and ZAMS stand for turn-off age, half-age, and zero-age main- 
sequence respectively. In the fifth column we show the number 
of SPH particles in the simulations; in all cases, we use well over 
100k equal mass SPH particles in order to achieve high resolution 
in the dense stellar interiors and to insure convergent results. 



stellar evolution code (|Paxtonll2004l '). In particular, we have 
used EZ to evolve a non-rotating zero age main-sequence 
(ZAMS) model with a primordial helium abundance Y ~ 0.3 
and metallicity Z — 0.02. The resulting star is composed of 
spherical mass-shells which store necessary stellar data, such 
as local composition, temperature, and density, as well as 
enclosed stellar mass and radius. The stars are evolved until 
the primary star reaches its half-age of its main-sequence 
life-time (HAMS) or until the primary has only 1% hydrogen 
abundance in its core (TAMS). In all cases the secondary 
star is evolved to the same age as the primary with which it 
collides. 



3.2 Results 

We have shown in Sect. [5]that the two main quantities which 
determine the structure of the collision product is buoyancy 
and composition. Since in main-sequence stars the mean 
molecular mass is mostly dependent on hydrogen mass frac- 
tion, we present in Fig. [T]buoyancy and hydrogen mass frac- 
tion profiles for a selection of the parent stars. In the case 
of TAMS stars, the primary consumed nearly all of hydro- 
gen in its core. Nevertheless, the secondary has consumed 
less than the half of its initial hydrogen supply. HAMS par- 
ents, on the other hand, are barely evolved; all of the parent 
models have burnt less than the half of the initial amount 
of hydrogen in their cores. 

The buoyancy of the parent stars is shown in the two 
lowest panels of Fig. [T] Even though the 2OM and 40M Q 
stars are at the end of the main-sequence, the buoyancy in 
their cores is noticeably higher than that of the secondary 
star. Therefore, based on the results of Sect. [2] it is natural 
to expect that the 8M© star will occupy the centre of the 
merger product. In the case of the IOMq primary, however, 
one may expect mixing, since the buoyancy of the primary 
and the secondary stars differ by less than a factor of two, 
and this difference can be modified by shock heating. 

In this work, we model shock heating by increasing the 
buoyancy of each fluid element. For this we need to assume 
a particul ar distribution of shock heating. Motivated by the 
results of lLombardi et al.l (|2002l ') , we adopt the following 
simple distribution: 



log(A f /Ai) = a + b ■ log(Pi/P c ,i) 



(10) 
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Figure 1. Internal structure of a selection of TAMS (left side) 
and HAMS (right side) parent stars which were used in simula- 
tions. The upper panels show density profile as a function of en- 
closed mass. The black solid and dotted lines on the left side show 
density profile of the 8Mq star from T18 and T48 models respec- 
tively, while on the right side the 8Mq star is taken from H18 
and H48 models. The red dash-dotted, green dash-dot-dot-dot 
and blue dashed lines show density profiles of IOMq , 20Mq and 
40Mq stars respectively. The second panels from the top show 
temperature profiles, while the third and fourth panels present 
hydrogen mass fraction and buoyancy profiles respectively. 

Here, Ai and Af are the initial and the final buoyancy of a 
fluid element, respectively, Pi is the initial pressure, and P c .i 
is the central pressure of a parent stars from which the fluid 
has come. We compute parameters a and b by fitting Eq. \TU\ 
to the simulation data. 

One of the advantages of using Eq. [10] can be seen by 
looking at Eq. [3] As an SPH particle preserves its compo- 
sition in the course of a simulation, Eq. [TO] describes the 
change in entropy due to the shock heating process. This, 
in turn, potentially expands the range of applicability of 
the sorting algorithm presented in Sect. [5] to an arbitrary 
equation of state, such as of degenerate matter inside white 
dwarfs. 

For every simulation from Tab. [JJ we obtain a and b 
values in Eq. [TUJ for each of the parent stars. We use this 
equation to compute the change in the buoyancy for each of 
the mass shells of the parent stars, and then we apply the 
sorting algorithm to generate the structure of the collision 
product which we compare with the simulation results. 

In Fig. [2] we present the results of the sorting algo- 
rithm and simulations. The sorting algorithm is able to re- 
produce two distinct branches in temperature, composition 
and buoyancy for TAMS collisions; the density, on the other 
hand, is single valued as one would naturally expect. In the 
case of HAMS collisions, the branches are less pronounced, 
as expected, because the parent stars have roughly the same 
composition (see Fig. [TJ. A noticeable discrepancy in the 
density of the outer layers is caused by an inadequate rep- 
resentation of shock heating in these layers by Eq. 1101 in- 
deed, both the density and buoyancy of our simple model 
diverge from the SPH results at about the same enclosed 



mass. The mass loss in the collisions from Tab. [TJ never ex- 
ceeds 10%: T18(8.3%), T28(8.9%), T48(5.0%), T88(1.9%), 
H18(6.8%), H28(4.7%), H48(2.1%), H88(0.8%), Z18(6.4%) 
and Z28(4.5%). As we expect, mass loss is less in collisions 
involving stars of significantly different mass, as less kinetic 
energy is redirected into the ejecta. In addition, we observe 
that mass loss is larger in collisions with evolved stars, due 
to their weakly bound envelopes. 

We note that the very existence of two branches in the 
results of Fig.[2]means that the buoyancy A does not increase 
outward at all locations. In the case of model T18, it is actu- 
ally the higher buoyancy branch, corresponding to hydrogen 
depleted fluid from the TAMS primary, which extends into 
the core of the merger product. Such a configuration is sta- 
ble only because it is accompanied by a sufficiently steep 
gradient in composition. Indeed, in model H18, the parent 
stars are too young to yield large gradients in composition, 
and consequently it is the lower buoyancy fluid from the sec- 
ondary that sinks to the centre of the merger product. It is 
also worth noting that in models Z18 and Z28, in which the 
ZAMS parent stars have uniform composition, we find that 
the buoyancy increases outwards throughout the collision 
products, as expected from Eq. [9] 

As we increase the resolution of parent stars as treated 
by the sorting algorithm, adjacent mass shells in different 
branches become increasingly closely spaced. The resulting 
large gradients will eventually be smoothed by the micro- 
scopic mixing between neighbouring fluid elements. To sim- 
ulate such mixing, we process the results from the sorting 
algorithm into equal-mass bins, such that each bin accom- 
modates a sufficiently large number of mass shells. We then 
average the stellar data, such that the total volume and ther- 
mal energy of each bin, as well as the total mass of each of 
the chemical elements, are conserved. From the average val- 
ues of density, thermal energy and composition, we compute 
the rest of the thermodynamic quantities, such as temper- 
ature. We find that this procedure converges once the high 
resolution model has roughly 10 5 mass shells averaged into 
200 bins. The simulation data is averaged in the same way. 

In Fig. |3] we show the results of the mixing applied to 
the models shown in Fig. [2] Whereas no significant mixing 
occurs in models T48 and H48 as the 8Mq star settled in 
the centre of the primary, the interior of the collision prod- 
uct in model T18 is mixed. In model H18, however, only 
the core of the secondary star settles in the centre of the 
primary star; the rest of the fluid is mixed. The resulting 
buoyancy increases outwards throughout the product in all 
models, except in the very centre of T18. Here, the decrease 
in buoyancy at ~ O.8M0 is accompanied by an increase in 
mean molecular mass in such a way that Eq. [9] is satisfied. 



4 DISCUSSION AND CONCLUSIONS 

We have presented a new sorting algorithm which generates 
the structure of collision products in hydrostatic equilib- 
rium. We have tested the algorithm by carrying out a set of 
SPH simulations of head-on stellar collisions between mas- 
sive main-sequence stars of different masses and ages. By 
calibrating the shock heating with the simulation data, we 
have been able to quickly generate the structure of the colli- 
sion products consistent with the results of the simulations. 
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Figure 2. The non-mixed structure of the collision products from a selection of simulations. The left set of plots displays results for 
TAMS collisions, whereas the right set displays the same data for HAMS collisions. The black dots are the results of simulations, while 
the green circles are the results of the sorting algorithm. The upper panels show the density profile as a function of enclosed mass for 
T18, T48, H18 and H48 simulation models from left to right respectively. In the second panels from the top, we show the temperature 
profile, and in the third and the fourth panels, we show the hydrogen mass fraction and buoyancy profiles respectively. 
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Figure 3. The structure of the mixed collision products from 
the same selection of models as presented in Fig. [2] Here, the 
black solid and red dotted lines show the results of the sorting 
algorithm and of the simulation respectively for the cases T18 on 
the left side and H18 on the right side. The blue dash-dotted and 
green dashed lines show the results of sorting algorithm and of 
the simulation respectively for the cases T48 on the left side and 
H48 on the right side. The upper panels show density profiles 
as a function of enclosed mass, and the second panel from the 
top shows temperature profiles. The third and the fourth panels 
from the top show hydrogen mass fraction and buoyancy profiles 
respectively. 



It takes less than a minute to generate a collision prod- 
uct for any of the models presented in Tab. [T] whereas an 
SPH simulations takes at least a day of CPU time on the 
same computer. Such a speed-up will help to make it possi- 



ble to include realistic stellar collisions in simulations of star 
clusters. 

We have found that the chemical composition and tem- 
perature can be a multivalued function of the enclosed mass, 
if microscopic mixing processes are ignored. This is exhibited 
by the fact that the neighbouring mass shells have a discon- 
tinuously changing temperature and chemical composition, 
whereas the density and pressure remains continuous. Us- 
ing the sorting algorithm, we have found that the situation 
is unaffected by increasing the number of shells. Instead, 
the chemical composition and temperature become double- 
valued in the continuous limit. Such an obviously unphysical 
situation is mended by microscopic mixing. We have simu- 
lated the mixing process by conservatively averaging multi- 
ple mass shells into a single mass bin. The resulted profiles 
are both single valued and consistent with the simulation 
data. 

We have found significant mixing in only the T18 and 
H18 models. The reason for this is that buoyancy in the 
core of the primary and the secondary stars are similar, 
whereas for the other models, such as T48 and H48, the 
buoyancy of the primary stars is substantially larger than 
that of the secondary stars. Therefore, in those cases, fluid 
from the secondary star occupies the central region of the 
collision product without being mixed with the fluid from 
the primary star. We therefore conclude, that hydrodynamic 
mixing in collisions between main-sequence stars occurs only 
when the masses and ages of the parents stars are similar to 
each other. We find that mixing already becomes unimpor- 
tant at the mass ratio q ~ 0.4, which is the mass ratio of 
models T28, H28 and Z28. In such cases, the secondary star 
simply occupi es the centre of the collision product wi thout 
much mixing ijDale fc Davies 2006; S uzuki et al.ll2007l ). 

Microscopic mixing may play an important role dur- 
ing the collision process itself. If the diffusion processes are 
ignored, the turbulent motion of the fluid can result in a 
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discontinuous distribution of chemical composition and tem- 
perature, while maintaining continuity in density and pres- 
sure. Such a distribution is susceptible to microscopic mixin g 
which tends to smooth the discontinuities (jDimot akis 2005). 
If the timescale of microscopic mixing is larger than the 
time a collision product takes to settle into quasi-hydrostatic 
equilibrium, the mixing can be ignored during the collision 
event and applied only as a post-collision process. On the 
other hand, if the timescale of the microscopic mixing is 
shorter than the timescale of a collision event, then one must 
include the effects of mixing into the simulation itself. 

The sorting algorithm can in principle be applied to 
any type of stars; for example, it is possible to apply it to 
the merger of white dwarfs. The basic requirement is simply 
to have the equation of state which gives the relationship 
between pressure, density, composition, and specific entropy 
(or buoyancy), analogous to Eqs. [3] and [4] While it is also 
necessary to determine mass loss and shock heating of the 
fluid, these can be estimated through energy conservation 
arguments. 
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